Three-Dimensional  Rayleigh-Benard  Convection  of  a 
Rarefied  Gas:  DSMC  and  Navier-Stokes  Calculations 

S.  K.  Stefanov*,  V.  M.  Roussinov*,  C.  Cercignani 


Institute  of  Mechanics,  BASci.,  Acad.  G.  Bonchev  St)'.,  Block.  4,  1113  Sofia  Bulgaria 
Dipartimento  di  Matematica,  Politecnico  di  Milano,  p.zza  L.  da  Vinchi,  32,  20133  Milan,  Italy 


Abstract.  The  three-dimensional  Rayleigh-Benard  convection  of  a  rarefied  gas  is  studied  numerically  by  using  DSMC 
and  Navier-Stokes  finite  difference  methods.  We  present  the  results  obtained  in  a  three-dimensional  box  with  an  aspect 
ratio  A=Lx:Ly:Lz  at  a  fixed  wall  temperature  ratio  r=Tc/Th=0.1.  In  general,  the  three-dimensional  results  confirm  the 
validity  of  the  zone  of  convection  obtained  in  the  previous  author’s  two-dimensional  considerations  [1,2],  A  detailed 
analysis  has  been  performed  for  Knudsen  number  Rn=0.02  and  A=2:2:l,  A=3:3:l,  and  A=  4:4:1  and  different  Froude 
numbers  ranged  within  the  zone  of  convection.  As  result,  a  variety  of  different  stable  convection  patterns  in  form  of  rolls 
and  squares  has  been  obtained.  For  the  lowest  Knudsen  numbers  computed,  Kn=0.001  and  Kn=0.002  ,  Froude  number 
Fr=  1.0,  and  a  fixed  aspect  ratio  A=2:2:l,  an  irregular  (chaotic)  flow  in  form  of  unstable  polygonal  patterns  has  been 
observed. 


INTRODUCTION 

The  formation  and  development  of  convection  flows  in  a  fluid  confined  between  two  horizontal  parallel  plates 
with  the  bottom  plate  heated  from  below  is  a  classical  problem  known  as  the  Rayleigh-Benard  problem.  In  our 
previous  papers  [1,2]  we  have  studied  the  formation  of  convection  patterns  in  the  Rayleigh-Benard  (RB)  flow  of  a 
rarefied  gas  in  a  rectangular  two-dimensional  domain  by  using  two  different  numerical  approaches:  particle 
simulation  by  the  DSMC  method  and  Navier-Stokes  finite  difference  (FD)  calculations.  The  calculations  have  been 
performed  within  a  range  of  Knudsen  number  Kn  =  0.001-0.04  and  Froude  number  Fr=0.5-1.5xl03  and  a  fixed 
temperature  ratio  r=0.1.  We  have  observed  a  variety  of  possible  final  states  (attractors),  namely:  pure  heat 
conduction,  stable  vortex  roll  convection,  periodic,  wavy,  and  chaotic  regimes. 

In  the  present  paper  we  extend  the  study  of  the  problem  to  convection  gas  flows  in  a  three-dimensional  box.  The 
low  Knudsen  numbers  allow  again  the  problem  to  be  investigated  by  using  two  numerical  approaches:  DSMC 
method  (particle  approach)  and  finite  difference  method  (continuum  approach  based  on  the  model  of  compressible 
viscous  heat  conducting  gas  with  state-dependent  transport  coefficients).  Using  the  DSMC  method,  the  three- 
dimensional  convection  of  a  rarefied  gas  has  been  investigated  by  Watanabe  [3]  for  the  case  of  a  hypothetical  value 
of  the  gravitational  acceleration  (the  Froude  number  in  non-dimensional  formulation)  so  as  to  minimize  the  density 
variation  in  the  pure  heat  conduction  state.  The  same  approach,  imitating  the  conditions  of  the  Boussinesq 
approximation  (valid  for  continuum  fluid  dynamics),  has  been  used  previously  by  Garcia  et  al.  [4]  and  Watanabe  et 
al.  [5]  for  the  two-dimensional  case  of  convection  in  rarefied  gas.  In  our  calculations,  the  governing  non- 
dimensional  parameter  Fr  was  varied  freely  so  strongly  stratified  states  could  be  investigated.  This  idea  was 
proposed  by  Stefanov  and  Cercignani  [6]  and  used  by  Sone  et.  al.  [7],  Stefanov  et.  al.  [1,2]  for  a  detailed 
investigation  of  the  two-dimensional  Rayleygh-Benard  problem  of  a  rarefied  gas.  In  the  present  paper  we  follow  that 
second  approach  and  investigate  the  three-dimensional  RB  flow  for  a  set  of  Kn  and  Fr  at  fixed  temperature  ratio 
r=0.1.  The  dependence  of  the  pattern  formation  on  the  aspect  ratio  A  =  Lx:Ly:Lz  is  also  studied.  The  most  detailed 
analysis  is  performed  for  Kn=0.02,  Fr=3.0  and  aspect  ratios  A=2:2:l,  A=3:3:l  and  A=  4:4:1.  For  the  same  Kn=0.02 
and  aspect  ratios  A=2.0  and  A=3.0  the  calculation  are  performed  for  Fr=1.5,  Fr=3.0  and  Fr=7.0.  The  next 
computations  have  been  performed  for  sufficiently  low  Knudsen  numbers  Kn=0.002  and  Kn=0.001,  where  one  can 
expect  development  of  more  complicated  flow  regimes,  as  was  discovered  in  [1,2]  for  the  two-dimensional  RB 
flow. 
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In  general,  the  analysis  of  the  obtained  results  for  the  three-dimensional  Rayleigh-Benard  flow  show  that  the 
zone  of  governing  parameters,  in  which  a  transition  from  pure  heat  conduction  to  the  regime  of  convection  takes 
place,  is  in  good  agreement  with  the  zone  delineated  by  the  two-dimensional  calculations  [1].  A  chaotic  convection 
in  form  of  different  unstable  polygonal  patterns  has  been  obtained  within  the  range  of  Knudsen  numbers  Kn=0.001- 
0.002  and  a  Froude  number  around  Fr=l. 0-1.1.  For  larger  Kn  and  Fr  we  have  observed  a  variety  of  stable  structures 
in  form  of  roll  or  square  patterns. 

Problem  Formulation  And  Computational  Considerations 

The  three-dimensional  Rayleigh-Benard  problem  for  a  rarefied  gas  treats  a  monatomic  simple  gas  with  average 
number  density  n0  studied  in  a  rectangular  computational  domain  D(LxxLyxLz)  i.e.  with  an  aspect  ratio 
A=Lx/Lz:Ly/Lz:Lz/Lz.  The  domain  is  confined  in  the  vertical  z-direction  by  two  diffusely  reflecting  horizontal  walls 
(the  lower  wall  temperature  is  greater  than  the  upper  one,  Th  >  Tc)  and  a  periodic  boundary  condition  is  implemented 
in  the  horizontal  x-  and  y-directions  at  vertical  planes  at  mutual  distances  Lx  and  Ly.  A  constant  acceleration 
g=(gx,gy,gz)=(0,0,-g)  acts  on  the  gas  in  each  point  of  the  computational  domain.  A  hard  sphere  molecular  model  is 
employed  in  the  simulation  of  binary  collisions.  The  governing  nondimensional  parameters  are  the  Knudsen  number 
Kn=f(/Lz,  the  Froude  number  Fr=Vth2/gLz,  and  the  temperature  ratio  r=Th/Tc,  where  €0  is  mean  free  path  in  a  “hard 
sphere”  gas  with  density  n0,  and  Vth  is  the  most  probable  molecular  velocity  in  equilibrium  gas  with  temperature  Th. 
In  our  consideration  all  variables  are  normalized  by  using  the  following  scales:  for  density,  p0=mn0  (m  is  the 
molecular  mass);  for  velocity,  Vth;  for  length,  the  distance  Lz;  for  time,  t0=Lz/Vlh;  for  temperature,  Th. 

The  time  evolution  of  this  system  when  starting  from  a  certain  initial  state  has  been  studied  by  means  of  the  two 
approaches  discussed  in  the  introduction:  molecular  and  continuum. 

The  DSMC  technique  uses  a  finite  number  of  model  particles  that  move  and  collide  in  the  three-dimensional 
space  domain  D  to  perform  a  direct  simulation  of  the  molecular  gas  dynamics.  The  algorithm  follows  the  basic  steps 
of  the  "No  Time  Counter"  scheme  suggested  by  Bird  [8].  A  uniform  rectangular  basic  grid  with  NxxNyxNz  cells 
covers  the  computational  domain.  Within  every  time  step  each  basic  cell  is  subdivided  dynamically  into  subcells  in 
order  to  meet  the  special  resolution  requirements;  the  size  of  a  subcell  is  less  than  the  local  mean  free  path.  In  this 
way,  only  the  information  about  the  basic  grid  is  kept  permanently  in  the  computer  memory.  The  sampling  is 
conducted  on  a  coarser  grid  with  larger  cells  containing  a  group  of  basic  cells.  This  multilevel  grid  scheme  allows  a 
correct  calculation  of  the  molecular  process  on  an  adaptive  fine  grid  and  at  the  same  time  provides  a  meaningful 
sample  size  of  the  instantaneous  macroscopic  variables  by  averaging  simultaneously  over  a  group  of  cells.  For 
example,  for  Kn=0.02  and  A=2:2:l  a  grid  with  100x100x50  basic  cells  is  used;  while  for  Kn=0.02  and  A=4:4:l  — 
200x200x50.  For  the  computed  lower  Kn=0.002  the  grid  at  A=2:2:l  is  with  200x200x100  basic  cells.  The 
instantaneous  fields  are  sampled  by  time  averaging  over  500  time  steps  At  (At  is  less  than  mean  collision  time) 
Finally,  the  obtained  three-dimensional  instantaneous  macroscopic  fields  are  processed  by  a  simple  average  filter  in 
order  to  remove  the  high-frequency  statistical  fluctuations  from  the  simulation  results.  Depending  on  the  Knudsen 
number  and  the  aspect  ratio  A  the  total  number  of  model  particles  varies  from  5xl06  to  32xl06. 

The  3D  finite  difference  (FD)  calculations  of  the  continuum  model  based  on  the  Navier-Stokes  equations  for  a 
compressible  viscous  gas  with  state  dependent  transport  coefficients  (see  [1])  are  performed  by  using  an  implicit 
conservative  finite  difference  (second  order  in  space  and  first  order  in  time)  scheme  with  inner  iteration  process 
within  each  time  step.  During  each  iteration  a  7-bandmatrix  linear  equation  system  obtained  for  each  macroscopic 
variable  is  solved  by  using  a  generalized  conjugate  gradient  (GCG)  method.  The  iteration  error  is  chosen  to  be  not 
larger  than  10"10  —  10"11  depending  on  the  chosen  computational  parameters  of  the  scheme.  The  size  of  the 
rectangular  uniform  computational  grid  depends  on  the  Knudsen  number  and  the  aspect  ratio  A.  For  the  most 
difficult  cases,  Kn=0.002  and  Kn=0.001  at  fixed  aspect  ratio  A=2.0,  a  grid  with  200x200x100  nodes  is  used.  For  the 
larger  Kn=0.02  the  calculations  are  performed  on  grids  100x100x50  ,  150x150x50,  and  200x200x50  for  aspect  ratios 
A=2.0,  A=3.0,  and  A=4.0,  respectively. 


Numerical  Results 

We  illustrate  the  obtained  results  using  a  limited  set  of  the  computational  data.  In  our  previous  works  [1,2]  we 
have  shown  that  the  two-dimensional  calculations  of  the  RB  flow,  performed  by  using  both  molecular  and 
continuum  approach  for  the  range  of  the  governing  parameters  mentioned  in  the  introduction,  exhibited  a  qualitative 
and  quantitative  similarity  of  the  obtained  convective  flows.  That  is  why,  here  our  main  goal  is  not  to  present  a 
detailed  comparison  of  both  model  calculations  for  the  three-dimensional  case  of  RB  flow  but  to  use  each  of  the  two 
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models  where  it  is  more  effective  in  order  to  cover  a  wider  range  of  the  parameters  where  new  interesting  convection 
phenomena  might  exist.  A  rough  preliminary  comparative  analysis  has  convinced  us  that  we  could  not  expect  great 
surprises  and  both  models  should  give  similar  results  within  the  studied  range  of  the  parameters.  This  work  is  in 
progress  and  a  more  detailed  discussion  on  the  comparison  of  the  two  models  in  the  three-dimensional  case  of  the 
RB  convection  will  be  given  in  a  separate  paper.  Here  we  follow  the  line  of  presenting  our  results  for  a  set  of 
calculations  performed  by  using  the  DSMC  simulations  for  Kn=0.02  and  various  Fr  and  aspect  ratios  where  the 
particle  method  was  more  effective  than  the  finite  difference  computations.  In  the  cases  with  a  smaller  Knudsen 
number,  Kn=0.002  and  0.001,  the  FD  method  was  considerably  more  effective  because  it  turned  out  possible  to 
obtain  a  quantitatively  correct  results  on  sufficiently  coarser  grid  than  a  grid  that  meets  the  computational 
requirements  of  the  DSMC  method.  Thus,  the  results  presented  here  for  the  chaotic  convection  at  Kn=0.002  and 
Fr=1.0  were  computed  by  using  the  FD  method.  This  case  was  run  by  DSMC  within  the  initial  part  of  the  transient 
period  where  the  results  were  similar  to  those  obtained  by  FD.  For  the  comparative  analysis  in  this  case,  the  DSMC 
calculations  will  be  continued  beyond  the  transient  period. 


FIGURE  1.  The  zone  of  convection  and  the  points  of  the  three-dimensional  computations:  the  analytical  conditions  for  the  left 
and  right  boundaries  of  the  zone  of  convection  [l].are  given  with  solid  lines,  the  gray  shaded  zone  illustrates  the  two-dimensional 
calculations  [1.2],  the  three-dimensional  calculations  are  given  with  filled  circles. 

Figure  1  shows  the  points  (filled  circles)  of  the  three-dimensional  calculations  in  a  coordinate  system  of  Fr  and 
Kn  numbers  (the  temperature  ratio  is  fixed  r=0.1).  The  gray  shaded  polygon  illustrates  the  zone  of  convection 
obtained  from  the  two-dimensional  calculations  [1],  The  analytical  conditions  for  the  left  and  right  boundaries  of  the 
zone  of  convection  (  eqs.  (15)  and  (16)  in  [1])  are  shown  by  solid  lines.  Additionally,  several  3D  calculations  have 
been  performed  for  points  outside  of  the  delineated  zone  of  convection  but  close  enough  to  its  boundaries.  In 
general,  these  3D-runs,  performed  by  using  DSMC  method  for  Kn=0.02  and  Kn=0.01,  confirmed  the  two- 
dimensional  results  [  1  ]  obtained  for  the  boundaries  of  the  zone  of  convection. 

Figures  2  and  3  illustrate  the  convection  patterns  for  fixed  Kn=0.02  and  A=2.0  at  Fr=1.5  (a)  (  the  left  side 


(a)  (b) 

FIGURE  2.  Velocity  vector  field  in  xy-plane  at  z=0.9for  Kn=0.02,  A=2,  and  Fr=1.5  (a)  ,  Fr=7.0  (b). 
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of  the  convection  zone  where  gravity  prevails  against  the  temperature  gradient  action)  and  at  Fr=7.0  (  the  right  side, 
where  the  temperature  gradient  is  dominating).  The  “view  form  above”  at  z=0.9  (Fig.  2)  showed  that  in  both  cases 
the  convection  flow  is  organized  as  an  asymmetrically  positioned  single  square  cell  with  a  single  sink  point  (the 
focal  point  of  the  velocity  vectors).  Figure  3  shows  the  similar  vortex  configuration  in  a  vertical  plane  section  at 
x=0.5  in  cases  (a)  and  (b)  as  defined  in  Figure  2.  From  the  figure  one  can  see  that  the  vortex  center  is  shifted  up  with 
the  increase  of  the  Froude  number. 


(a)  (b) 

FIGURE  3.  Velocity  vector  field  in  yz-plane  at  x=0.5for  Fr=1.5  (a)  and  Fr  =7.0  (b) 

Figures  4  and  5  show  the  flow  patterns  for  different  aspect  ratios  A=2,  A=3,  and  A=4  at  Kn=0.02  and  Fr=3.0. 


y 


FIGURE  4.  Velocity  vector  field  in  xy-plane  at  z=0.9  for  Kn=0.02  ,  Fr=3.0  and  aspect  ratios  A=2  (left),  A=3  (right) 


FIGURE  5.  Velocity  vector  field  in  xy-plane  at  z=0.9  (left)  and  3D  velocity  vector  plot  for  Kn=0.02  ,  Fr=3.0, 

A=4.0. 
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The  flow  pattern  for  A=2  is  again  a  single  square  cell.  For  A=3  the  result  is  different  —  inclined  square  cells  with 
a  size  equal  to  the  pattern  size  in  the  case  A=2.  The  picture  for  A=4  is  drastically  changed:  the  flow  patterns  are 
stable  four  rolls  (Fig.  5,  (right))  with  hardly  seen  in  Figure  5  (left)  variable  thickness  chess  positioned.  This  effect 
can  be  seen  from  the  density  plot  at  the  middle  plane  z=0.5  in  Fig.  6  (left).  In  order  to  keep  this  configuration  stable 
the  velocity  field  is  structured  in  a  complicated  way  illustrated  in  Fig. (6)  (right)  with  the  xy-velocity  field  at  z=0.5. 
Here  one  can  observe  a  stable  low  intensity  secondary  flow  in  xy-  directions. 


FIGURE  6.  The  density  field  at  z=0.5  (left)  and  the  3D  velocity  vector  plot  (right)  for  Kn=0.02,  A=4. 

The  analysis  of  the  flow  patterns  obtained  by  DSMC  simulations  for  Kn=0.02  leads  to  following  two 
conclusions.  First,  the  flow  configuration  depends  strongly  on  the  boundary  conditions  and  the  aspect  ratio.  Thus, 
for  the  small  aspect  ratio  A=2  we  have  a  single  square  cell  for  all  magnitudes  of  Fr  ranged  within  the  zone  of 
convection.  For  A=3  the  patterns  are  the  same  in  form  and  size  but  are  inclined.  Similar  inclined  patterns  have 
obtained  by  Watanabe  [3]  for  odd  aspect  ratios.  This  can  be  explained  with  the  link  between  the  periodicity 
boundary  conditions  and  the  odd  aspect  ratio.  Second,  there  are  long  ranged  correlations  which  are  cut-off  in  the 
small  aspect  ratio  computations.  This  conclusion  is  supported  by  the  results  shown  in  figures  5  and  6.  The  case  with 
A=4  shows  that  the  larger  computational  box  allows  four  roll  configuration;  this  is  not  observed  in  the  cases  with 
A=2.0  and  A=3.0. 

Figures  7  (a)  —  (f)  present  a  case  of  chaotic  Rayleigh-Benard  convection  at  Kn=0.002,  Fr=1.0,  and  r=0.1, 
obtained  by  using  Navier-Stokes  calculations.  The  series  of  snapshots  taken  at  different  times  after  the  end  of  the 
transient  period  illustrates  the  time  evolution  of  the  convection  patterns  within  a  period  of  about  100  units  in  the 
chosen  time  scale.  The  times  are  selected  in  a  way  that  the  xy-velocity  vector  fields  at  a  fixed  z=0.9  (close  to  the 
cold  wall)  show  a  variety  of  typical  flow  patterns.  Most  of  them  are  in  the  form  of  unstable  polygonal  configurations 
interacting  with  each  other  in  a  complicated  way.  In  the  figures  the  narrow  stripes  with  gathering  arrows  of  the 
velocity  vectors  delineate  the  edges  of  the  polygons  where  the  cold  flux  goes  down  from  the  cold  (top)  to  the  hot 
(bottom)  wall.  Within  a  wide  area  inside  the  polygons  the  flux  is  directed  up.  The  size  of  the  configurations  varies 
from  small  flow  patterns  (Fig.  7  (c))  less  than  1/10  of  the  computational  box  size  to  large  ones  (Fig.  7  (e)) 
comparable  to  it. 
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(e) 

FIGURE  7.  Instantaneous  velocity  vector  filed  for  Kn=0.002,  Fr=l, 
t=220.0;  (c)  t=235.0;  (d)  t=260.0;  (e)  t=280.0;  (f)  t=  300.0. 


(f) 

A=2  at  different  times:  (a)  t=200.0;  (b) 
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